Kluczowym elementem teorii produkcji energii termojądrowej w gwiazdach i nukleosyntezy jest umiejętność obliczenia lub zmierzenia tempa zachodzenia reakcji. W przypadku ,,chmury’’ cząstek o gęstości  zderzającej się z ,,chmurą’’ cząstek o gęstości liczbowej  poruszającą się z prędkością względną tempo zachodzenia reakcji   (liczba zderzeń na jednostkę czasu w jednostce objętości) wynosi:



gdzie  to przekrój czynny na zderzenie. Gdyby zderzały się cząstki sztywne, np. kule, przekrój czynny byłby po prostu sumą geometrycznych przekrojów czołowych. W przypadku cząstek kwantowych,  wynika z odpowiednich teorii fundamentalnych lub ich przybliżonych modeli kalibrowanych eksperymentami laboratoryjnymi.

W przypadku gwiazd, reagujące składniki poruszają się losowo zgodnie z prawami fizyki statystycznej. Dla Słońca i większości gwiazd wystarczającym przybliżeniem jest rozkład Maxwella. Aby ułatwić Mathematice całkowanie, oznaczamy jednym symbolem kT energię termiczną .

$Assumptions={kT>0,m>0}

Poniżej mamy standardowy rozkład Maxwella jednej ze składowych prędkości częstek gazu.

Integrate[4Pi v^2 Exp[-m v^2/(2 kT)],{v,0,Infinity}]

Rozkład normalizujemy tak aby całka dawała gęstość liczbową cząstek .

n Sqrt[m/(2 Pi kT)]Integrate[ Exp[-m vx^2/(2 kT)],{vx,-Infinity,Infinity}]//Simplify

W przypadku zderzeń w gazie, wzór na tempo reakcji przyjmuje postać całki.

Powyższą całkę poszóstną najprościej obliczyć przechodząc do układu środka masy (CM). To ponadto ułatwia wyizolowanie , które zwykle jest mierzone lub oblicane właśnie w układzie CM.

sysCM={
vCMx==(m1 v1x+m2 v2x)/(m1+m2)
,
vCMy==(m1 v1y+m2 v2y)/(m1+m2)
,
vCMz==(m1 v1z+m2 v2z)/(m1+m2)
,
vRelx==v1x-v2x,vRely==v1y-v2y,vRelz==v1z-v2z};

varOLD={v1x,v2x,v1y,v2y,v1z,v2z};

zamianaZmiennych = Solve[
sysCM,varOLD
]

tr=Last/@sysCM

varNEW=First/@sysCM

Szczęśliwie dla wielu studentów, Jakobian transformacji wynosi 1.

Table[D[tr[[i]],varOLD[[j]]],{i,1,6},{j,1,6}]//Det//Simplify

Prosty test poprawności wzorów: energia w układzie CM.

1/2 m1 (v1x^2+v1y^2+v1z^2)+1/2 m2 (v2x^2+v2y^2+v2z^2)-1/2(m1+m2)(vCMx^2+vCMy^2+vCMz^2)/.zamianaZmiennych[[1]]

Simplify[%]

f[v_,m_]:=Sqrt[m/(2Pi kT)]Exp[-m v^2/(2 kT)] (* Definicja rozkładu Maxwella *)

i=Inactive[Integrate][
f[v1x,m1]f[v1y,m1]f[v1z,m1]f[v2x,m2]f[v2y,m2]f[v2z,m2] σ Sqrt[(v1x-v2x)^2+(v1y-v2y)^2+(v1z-v2z)^2],{v1x,-Infinity,Infinity},{v1y,-Infinity,Infinity},{v1z,-Infinity,Infinity},{v2x,-Infinity,Infinity},{v2y,-Infinity,Infinity},{v2z,-Infinity,Infinity}]

$Assumptions=kT>0&&m1>0&&m2>0&&(vCMx|vCMy|vCMz|vRelx|vRely|vRelz)∈(Reals)

IntegrateChangeVariables[i,{vCMx,vCMy,vCMz,vRelx,vRely,vRelz},{
vCMx==(m1 v1x+m2 v2x)/(m1+m2)
,
vCMy==(m1 v1y+m2 v2y)/(m1+m2)
,
vCMz==(m1 v1z+m2 v2z)/(m1+m2)
,
vRelx==v1x-v2x,vRely==v1y-v2y,vRelz==v1z-v2z}]

DeleteCases[%,{True,True}]

i2=(InactiveHead[Integrate])[(1/(8 (kT^3) (π^3)))(E^(-(((m1^2) ((vCMx^2)+(vCMy^2)+(vCMz^2))+(m2^2) ((vCMx^2)+(vCMy^2)+(vCMz^2))+m1 m2 (2 (vCMx^2)+2 (vCMy^2)+2 (vCMz^2)+(vRelx^2)+(vRely^2)+(vRelz^2)))/(2 kT (m1+m2))))) (((m1 m2))^(3/2)) (Sqrt[(vRelx^2)+(vRely^2)+(vRelz^2)]) σ ,{vCMx,-∞,∞},{vCMy,-∞,∞},{vCMz,-∞,∞},{vRelx,-∞,∞},{vRely,-∞,∞},{vRelz,-∞,∞}]

IntegrateChangeVariables[Inactive[Integrate][((m1*m2)^(3/2)*Sqrt[vRelx^2 + vRely^2 + vRelz^2]*σ)/
   (E^((m1^2*(vCMx^2 + vCMy^2 + vCMz^2) + m2^2*(vCMx^2 + vCMy^2 + vCMz^2) + m1*m2*(2*vCMx^2 + 2*vCMy^2 + 2*vCMz^2 + vRelx^2 + vRely^2 + vRelz^2))/(2*kT*(m1 + m2)))*
    (8*kT^3*Pi^3)), {vCMx, -Infinity, Infinity}, {vCMy, -Infinity, Infinity}, {vCMz, -Infinity, Infinity}],{V,Vθ,Vφ},"Cartesian"->"Spherical"]

integrand2=Activate[%]

IntegrateChangeVariables[Inactive[Integrate][integrand2,{vRelx, -Infinity, Infinity}, 
  {vRely, -Infinity, Infinity}, {vRelz, -Infinity, Infinity}],{V,Vθ,Vφ},"Cartesian"->"Spherical"]

Activate[%]

wynik=PowerExpand[%]

Block[{m1=EulerGamma,m2=Glaisher, kT=Catalan},Inactive[NIntegrate][
f[v1x,m1]f[v1y,m1]f[v1z,m1]f[v2x,m2]f[v2y,m2]f[v2z,m2] Sqrt[(v1x-v2x)^2+(v1y-v2y)^2+(v1z-v2z)^2],{v1x,-Infinity,Infinity},{v1y,-Infinity,Infinity},{v1z,-Infinity,Infinity},{v2x,-Infinity,Infinity},{v2y,-Infinity,Infinity},{v2z,-Infinity,Infinity},MinRecursion->12, MaxRecursion->18]
]

Activate[%]

Block[{m1=EulerGamma,m2=Glaisher, kT=Catalan},((2 (Sqrt[kT]) (Sqrt[m1+m2]) (Sqrt[(2/π)]) )/((Sqrt[m1]) (Sqrt[m2])))]//N

%/%%-1

μ= (1/m1+1/m2)^-1

((2 (Sqrt[kT]) (Sqrt[m1+m2]) (Sqrt[(2/π)]) )/((Sqrt[m1]) (Sqrt[m2])))/(2 Sqrt[2 kT/Pi/μ])

Simplify[%]

Sqrt[8/Pi/μ](kT)^(-3/2)Integrate[En Exp[-En/kT],{En,0,Infinity}]

%/wynik//Simplify

Inactive[Integrate][(((kT*m1*m2)/(m1 + m2))^(3/2)*Sqrt[2/Pi]*V^3*σ[1/2 μ V^2])/(E^((m1*m2*V^2)/(2*kT*m1 + 2*kT*m2))*kT^3), {V, 0, Infinity}]

IntegrateChangeVariables[%,En,En==1/2 μ V^2]

%/.m1 m2->μ (m1+m2)

Simplify/@%

%/.σ->Function[En,S[En]Exp[-2 Pi η]/En]